function [err]  = stdp_spike_dt(varargin)
global CONTROL
err = 0;
%try
    if(nargin > 0)
        sf = varargin{1}; % jjust use the calling argument to access the trace...
    else
        sf = getmainselection;
    end;
    if(sf > 0)
        pflag = getplotflag;
        QueMessage('epsp_spike_dt analysis', 1); % clear the que
        for i = 1:length(sf)
            err = EPSC_TC2(sf(i), pflag);
            if(err ~= 0)
                fprintf(1, 'Error stdp_spike_dt on file: %s%s\n', CONTROL(i).path, CONTROL(i).filename);
            end;
        end;
    end;
% catch
%     watchoff;
%     QueMessage('Error in epsp_spike_dt Analysis routine', 1);
% end
return;


function [err] = EPSC_TC2(sf, plot_flag)

global VOLTAGE CURRENT
global CONTROL DFILE
% first generate empty output arrays

err = 0;
EPSPresult=[];

QueMessage('epsc_spike_dt - Starting');
dat = [];
time = [];
ivresult = []; % initialize it
do_spike = 0;	% 0 turns off detailed spike analysis; 1 turns it on
depol = 0;

[DFILE, err] = analysis_setup(DFILE, sf);

if(err ~= 0)
    return;
end;

%-----------------------------Begin analysis--------------------------------

% abstract general information
spike_thresh=number_arg(CONTROL(sf).thresh); % to delete traces with spikes...
QueMessage(sprintf('epsp_spike_dt: Spike Threshold: %7.2f mV', spike_thresh));
protocol=deblank(lower(CONTROL(sf).protocol));
rate = DFILE.rate.*DFILE.nr_channel/1000;
[records,pts]=size(CURRENT);

%if(~isempty(CONTROL(sf).psp_reclist))
%	[RL, err] = record_parse(CONTROL(sf).psp_reclist);
%end;
[RL, err] = record_parse(CONTROL(sf).reclist);
frec=min(RL);
lrec=max(RL);

QueMessage('epsp_spike_dt - time base');
% compute the time base for plotting (time is [rec, npts] array for EACH record)
time=make_time(DFILE);
tmax=max(max(time));
k = find(DFILE.rate < 1);
DFILE.rate(k) = 50;
RATES = (DFILE.rate .* DFILE.nr_channel)/ 1000; % array of sampling rates, convert to msec


QueMessage('epsp_spike_dt - Smoothing');
for i = 1:records
    if(~isempty(DFILE.shkl)) % any shock latency information available?
        shkp = floor(DFILE.shkl/RATES(i));
        VOLTAGE(i,shkp:shkp+5) = VOLTAGE(i,shkp-1); % artifact suppression
    end;
    fsamp = 1000/RATES(i); % get sampling frequency
    fco = 1000;		% cutoff frequency in Hz
    wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
    if(wco < 1) % if wco is > 1 then this is not a filter!
        [b, a] = fir_win(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
        vsmo(i,:) = DigitalFilt(b, a, VOLTAGE(i,:)')'; % filter all the traces...
    else
        vsmo(i,:) = VOLTAGE(i,:);
    end
end

arate = mean(RATES);
trmp = floor(DFILE.durs(1)/arate); % for RMP/ihold determination
trmpb = floor((DFILE.durs(1)-1)/arate);
[trmpb trmp];
sl=[];
if(length(DFILE.durs) > 1)
    sl(1)=DFILE.durs(1);
    for i = 2:5
        sl(i)=sum(DFILE.durs(1:1+(i-1)*2));
    end;
    slb=sl-1;
    tsl = floor(sl/arate);
    tslb = floor(slb/arate);
    vdep = 0;
    for j = 1:records
        for i = 2:5
            vdep = vdep + mean(vsmo(j,tslb(i):tsl(i)));
        end;
        vrmp(j) = mean(vsmo(j, trmpb:trmp));
    end;
    depol = vdep/(records*4)-mean(vrmp);
else
    depol = mean(vsmo(1:records, trmpb:trmp));
end;


%vsmo = voltage; % in case you don't want filtering
QueMessage('epsp_spike_dt - finding spikes');

% find spikes in the records.
min_isi = 5;
[first_spike, first_isi, nr_spikes, spike_train] = find_spikes(DFILE, 0, 600, spike_thresh);
v = size(spike_train);
ISILAT=[]; % initialize the ISILAT structure
for i=1:v(2)
    ISILAT(i).nspikes = length(spike_train(i).latency);
    ISILAT(i).isi= diff(spike_train(i).latency);
    ISILAT(i).lat = spike_train(i).latency;
    ISILAT(i).source = spike_train(i).source; % save the source trace number too
end
if(isfield(CONTROL(sf), 'spike'))
    spike = CONTROL(sf).spike; % get currentn spike...
else
    spike = [];
end;
spikesumm = spike_sum(DFILE, vsmo, ISILAT, 0, CONTROL(sf).Rmp, CONTROL(sf).filename,...
    CONTROL(sf).recbeg, CONTROL(sf).deadwin, min_isi, 5, 10, 1);
% summary spike information:
if(~isempty(spikesumm))
    n = size(spikesumm.vmax);
    spike.timestamp = datestr(now);
    spike.edit = 0; % not editable...
    spike.n = n(2);
    spike.measthresh = mean(spikesumm.measthresh); % fraction of height from rest to spike peak where measurement was made
    athr = find(~isnan(spikesumm.vthr));
    if(isempty(athr))
        spike.vthresh=NaN;
        spike.vthrmin=NaN;
    else
        spike.vthresh = mean(spikesumm.vthr(athr)); % true threshold (2nd derivative method)
        spike.vthrmin = min(spikesumm.vthr(athr)); % min value (usually the first spike...)
    end;
    spike.deadwin = mean(spikesumm.deadwin);
    spike.minisi = mean(spikesumm.minisi);
    spike.max=mean(spikesumm.vmax); spike.max_sd = std(spikesumm.vmax);
    spike.apht = mean(spikesumm.apht); spike.apht_sd = std(spikesumm.apht);
    spike.half_ht = mean(spikesumm.half_ht); spike.half_ht_sd = std(spikesumm.half_ht);
    spike.min=mean(spikesumm.vmin); spike.min_sd = std(spikesumm.vmin);
    spike.dvmax = mean(spikesumm.dvmax); spike.dvmax_sd = std(spikesumm.dvmax);
    spike.dvmin = mean(spikesumm.dvmin); spike.dvmin_sd = std(spikesumm.dvmin);
    dvvmax = spikesumm.dvvmax(find(~isnan(spikesumm.dvvmax)));
    spike.dvvmax = mean(dvvmax); spike.dvvmax_sd = std(dvvmax);
    dvvmin = spikesumm.dvvmin(find(~isnan(spikesumm.dvvmin)));
    spike.dvvmin = mean(dvvmin); spike.dvvmin_sd = std(dvvmin);
    hwa = spikesumm.hwa(find(~isnan(spikesumm.hwa))); hwb = spikesumm.hwb(find(~isnan(spikesumm.hwb)));
    spike.hwa = mean(hwa); spike.hwa_sd = std(hwa); % eliminate the ones where coulnd't get a measurement
    spike.hwb = mean(hwb); spike.hwb_sd = std(hwb);
    spike.ahpdepth =  mean(spikesumm.ahpdepth); spike.ahpdepth_sd = std(spikesumm.ahpdepth);
    CONTROL(sf).spike=spike; % make it a substructure
    % spike train stuff...
    if(plot_flag >= 0)
        h = findobj('Tag', 'CC_IV_spiketrain'); % check for pre-existing window
        if(isempty(h)) % if none, make one
            h = figure('Tag', 'CC_IV_spiketrain', 'Name', 'spike train Analysis', 'NumberTitle', 'off');
        end
        figure(h); % otherwise, select it
        clf; % always clear the window...
        fsize = 7;
        msize = 3;
        subplot(2,2,1);
        plot(spikesumm.lat, spikesumm.dvmax, 'gx');
        xlabel('T (ms)');
        ylabel('dvmax (mV/ms) - rise');
        subplot(2,2,2);
        plot(spikesumm.lat, spikesumm.dvmin, 'rx');
        xlabel('T (ms)');
        ylabel('dvmin (mV/ms) - fall');

        subplot(2,2,3);
        plot(spikesumm.lat(find(~isnan(spikesumm.hwa))), hwa, 'go');
        xlabel('T (ms)');
        ylabel('hwa (ms) - rising half-width');
        subplot(2,2,4);
        plot(spikesumm.lat(find(~isnan(spikesumm.hwb))), hwb, 'ro');
        xlabel('T (ms)');
        ylabel('hwb (ms) - falling half-width');

    end;
end;

CONTROL(sf).spike=spike;

EPSPK.spkl = mean_var(first_spike);
EPSPK.spkn = sum(nr_spikes)/length(nr_spikes);
EPSPK.shkl = DFILE.shkl;
x=number_arg(CONTROL(sf).psp_time);
EPSPK.shkdt = EPSPK.spkl-(EPSPK.shkl+4); % -x(1); % measured latency! (must be done manually)
EPSPK.depol = depol; % store the "pre spike" depolarization, calculated according to Sjostrom et al.
EPSPK.freq = 1000/(DFILE.durs(2)+DFILE.durs(3)); % calculate the frequency of the shoek train...
CONTROL(sf).div = sprintf('%d', floor(EPSPK.freq));
CONTROL(sf).EPSPKdt=EPSPK;
%if(isempty(CONTROL(sf).Commentary) | isspace(CONTROL(sf).Commentary))
CONTROL(sf).Commentary=sprintf('dt = %7.2f', EPSPK.shkdt);
datac('Refresh');

QueMessage(sprintf('epsp_spike_dt: dt = %6.2f ms', EPSPK.shkdt), 1);
return;

